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Abstract 

The usefulness of Gabor frames depends on the easy computability of a suitable dual window. This question is addressed under 
several aspects: several versions of Schulz’s iterative algorithm for the approximation of the canonical dual window are analyzed 
for their numerical stability. For Gabor frames with totally positive windows or with exponential B-splines a direct algorithm 
yields a family of exact dual windows with compact support. It is shown that these dual windows converge exponentially fast to 
the canonical dual window. 


Introduction 

The discrete Gabor transform is a useful tool for the analysis and synthesis of nonstationary signals. It is based on the 
representation of the energy distribution of a signal in the time-frequency plane. Its applications range over the decomposition 
of musical and acoustical signals [1], wireless communication |[3l, El and to the analysis of EEG signals |[5], ISl- Eor a 
given window function g G L^(]R) and lattice parameters q:,/ 3 > 0, the system of all corresponding time-frequency shifts 

g{g, a, /3) = {Mip Tko. g = g{- - ka) \ kj G Z} 

is called a Gabor system for L^(]R). It is called a Gabor frame for L^(]R), if there exist constants A,B>0, such that 

A||/f < ^ K/,M,;3Tfc„p)p<i3||/f, V/gL2(M). (1) 

k,l& 

The constants A, B are called lower and upper frame bounds of G{g, a, P). If g{g, a, /3) fulfills only the right hand inequality, 
it is called a Bessel sequence and B a Bessel bound. It is known that the frame inequality O implies the existence of a dual 
Gabor frame g{'j, a, (3) with dual window 7 G L^(]R), such that every / G L^(]R) can be represented as 

/ = ^ (/, MipTk^g) MipTkcl- ( 2 ) 

Eor a given Gabor system, the Gabor ttansform of a signal / is defined as the analysis operator 

Cg ■. L^{^) ^ 

Cgf = ((/, Mii3Tkag))k,leZ- 

The coefficient (/, Mij^Tkag) represents the energy distribution of / near the point (fca, Z/3) in the time-frequency plane. It 
may also be interpreted as the amplitude of the frequency Z/3 at time ka, insofar as such an interpretation is compatible with 
the uncertainty principle. The associated synthesis operator for the reconsttuctions (|2]l is the adjoint operator C*, and the frame 
operator Sg : (M) —?> (K.) is defined by 

Sgf = c; Cgf = Y, if^ MipTke^g) MipTke.g. 
k,iez 

In general, there exist many dual windows suitable for the reconstruction (|2|i. The standard choice is the canonical dual 
window 7 ° = S~^g. Eor a characterization of all dual windows see Q, 0. Since the applicability and usefulness of Gabor 
frames depends heavily on the knowledge and computability of a dual window, the numerical construction of dual windows 
has motivated numerous studies. As representative contributions we mention 0, USl and the large time-frequency analysis 
toolbox (LTEAT) ifTTI . 

Our contribution to the analysis of dual Gabor windows is twofold. On a general level, we study numerically stable methods 
for the computation of the canonical dual window 7 ° = S~^g. On a specific level, we study the efficient construction and 
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the behavior of a sequence of dual windows for Gabor frames with totally positive window functions and with exponential 
B-splines. 

We first present two stable implementations of a conventional iterative algorithm to approximate 7 °. The algorithm was 
originally proposed by Schulz IfT^ for matrix inversion. It is based on the Neumann series for the inverse frame operator 
and converges quadratically, see Algorithm |2l We provide a detailed analysis of the numerical error and show that our two 
implementations (the operator version and the hrst vector version) are stable. By contrast, the implementation proposed by 
Janssen (see ca and m) is often unstable because the numerical error roughly doubles in each step. Therefore the first 
two implementations are much preferable. As two illustrative examples, we use the window functions MONSTER dehned in 
|9l and a Gaussian window g, in order to compare all three implementations. The numerical results agree precisely with the 
predicted behavior of the numerical error. 

We then describe recent results on special Gabor systems whose window function is a totally positive (TP) function of hnite 
type or an exponential B-spline (EB-spline). TP functions are remarkable because so far they are the only window functions 
for which a complete characterization of all lattice parameters such that Q{g,a,f3) is a frame is known. More precisely, the 
Gabor system Q{g,a,P), with a TP function g of hnite type > 2 is a frame if and only if < 1 US). Subsequently, 
similar arguments in am showed that the Gabor system Q{BA,a, [3) of an EB-spline constitutes a frame for a = 1, (3 < 1, 
and some other lattice parameters, too. The proofs also provide a constructive method for the computation of inhnitely many 
dual windows 71 , with compact support, which we summarize in Algorithm |5] for TP functions of hnite type and Algorithm [T] 
for EB-splines. 

This construction offers several new and useful aspects that are special for TP windows and not shared by general window 
functions. 

(i) Algorithms |5] and |7] provide a family of dual windows 71 , both in hnite and inhnite dimensional models, namely for 

continuous signals in L^(K.), for discrete signals in and for periodic discrete signals in C^. Currently available toolboxes, 

such as LTEAT HT], work only for hnite-dimensional signals. 

(ii) The dual windows 7 ^ possess compact support of size 0{L), whereas the canonical dual 7 ° is known to have inhnite 
support. 

(iii) The dual windows 7 /, are exact and satisfy (EJ. This is in contrast to the standard iterative methods for the approximation 
of the canonical dual (see e.g. Algorithm |2|i, which generate only approximations of a dual window. 

As our main mathematical result we prove that the dual windows 7 ^ are good approximations of the canonical dual window 
7 ° = Sg^g and we show that they converge exponentially fast to the canonical dual window, i.e., || 7 i — 7°||2 = 0{e~P^). 
Therefore, by specifying the parameter L, Algorithms |5] and [T] provide a dual window 71 , with compact support and which 
approximates the canonical dual at a desired rate. 

The proof uses some ideas of the non-symmetric finite section method, but also requires a new technique related to the 
formulation of the Moore-Penrose pseudo-inverse of inhnite matrices in terms of orthogonal projections. 

As our main numerical contribution, we study and implement the case of discrete Gabor frames. We present some fast and 
stable algorithms to evaluate and discretize TP functions and EB-splines and their dual windows computed by the Algorithms |5] 
and|2] These algorithms are proposed as extensions to the Large Time Erequency Analysis Toolbox described in ifTTll . 

The paper is organized as follows: In section J] we study the numerical stability of a fast iterative algorithm for the 
approximation of the canonical dual window. In section|II]we summarize the algorithms for the construction of dual windows of 
TP functions and EB splines. In section |III] we formulate and discuss the main theorems about the convergence of the compactly 
supported dual windows 7 z, to the canonical dual window 7 °. Section |IV] explains some details about the implementation of 
Gabor frames with TP functions and EB splines. The appendix contains the technical details of the proofs of the main results. 

I. Some iterative algorithms eor approximating the canonical dual 

In this section we describe two iterative algorithms for approximating the canonical dual of an arbitrary frame = {fj}j^i 
for a Hilbert space TL. The central part of such algorithms is the approximation of the inverse of the corresponding frame 
operator 

jei 

We discuss the convergence and the numerical stability of various implementations. Einally we present some numerical tests. 
The following approximation schemes are proposed in the literature. 

Algorithm 1 (Erame algorithm). Choose 0 < X < 2/B, with B the upper frame bound of J-. Then q := ||J— AS’jrH < 1 and 

CO 

n—0 
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The partial sums of this Neumann series can be computed iteratively by 

Kq = XI, 

Kk+i = XI + {I-XSr)Kk, fee No. (3) 

The convergence rate is of order —= 0{ci^). The k’th approximation of the canonical dual frame iF° = fj}jei 
is given by {Kk fj}jei- 

This algorithm is very robust, but slow if q is close to one. This algorithm can be accelerated ini with conjugate gradient 
techniques and with a convergence rate (after k iterations regardless of whether the frame bounds A, B are known 
or not. 

An even faster method goes back to Schulz ifT^ and Hotelling ifTSl . 

Algorithm 2 (Schulz iteration). Choose 0 < X < 2/B. The version of Schulz iteration with ’’initial scaling” / IT?) Algorithm IV] 
is 


Jo=XI, 

Jk+l = ‘2 Jk — Jk Jk, fc G Nq. (4) 

This iteration implies the identity Jk = K 2 k_i and is therefore connected to the frame algorithm. The Schulz algorithm 
converges quadratically, i.e. 

- Jk+l\\ < ll^^ll 11^^' - Jfcf = (5) 

This algorithm was first described by Schulz m who used this method for matrix inversion. 

Proof: The claims in Algorithm |2] are proved by induction. Since Schulz’s algorithm is not as known as other iterative 
algorithms, we sketch the main steps. We first show that 

I — Sj^Jk = I — JkSjr = (/ — . (6) 

Assuming that ® is correct for k G No, we obtain 

/ — SjrJk+l = I — Sj^{2Jk — JkSjrJk) 

= {I-SrJkf = ({I-XSrf)\ 

as claimed. Using we show again by induction that Jk = K 2 k_i = X'^jJ^ (/ — XSj^y: 

2fc+l_i 

A ^ (i-xs^y 

j=o 

2fc_i 2'=-! 

= A ^ - xs^y + (/ - xs^f A ^ (/ - xs^y 

j=0 j=0 

— Jk “t” (f JkSj^'jJk — Jk+1 . 

The quadratic convergence rate now follows from the convergence properties of the Neumann series. ■ 

We discuss the implementation of Algorithmin the case of a Gabor frame Q{g, a, /3). Recall that the canonical dual frame 
is determined by the dual window 7 ° = Sf^g. We compare three different implementations of the Schulz iteration and provide 
some heuristics for their numerical stability. 

(i) Operator form: The numerical computation of the Schulz iteration as stated in Algorithm|2]provides operators Jk = Jk + Ek, 
where Ek denotes the accumulated forward error. Let Yk+i denote the new roundoff error in the k + I’st iteration, then the 
operator after k + 1 iterations of (l4]l is 

Jk+l — 2 Jk Jk Jk T YkJ^-i 

= Jk + 1 + Ek {I — Sg Jk) + {I — Jk Sg) Ek + Yk+1 + Cl(||£^fc||^). 

Since I — Sg Jk = I — JkSg = {I — XSg)^'° , we have 

\\Ek+i\\ = \\Jk+i - Jk+i\\ < 2 / \\Ek\\ + lln+ill + 0{\\Ekf). 

This estimate shows that the error accumulated in the first k iterations is damped and only a new round-off error is added. 
Hence this iteration is numerically stable. 
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(ii) Vector form: We compute approximations of the dual window 7 ° directly by setting jk ■= Jk 9 with the following 
algorithm: 


70 = Ap 

7fc+i = 2 7fc - Cg 7fc, fc G No . (7) 


Proof: We use the fact that Sjr and thus all Jk commute with the time-frequency shifts MipTja. This commutation rule 
implies the identity 

^*k ^ ~ ^ Mip Tja g = JkCgC 

for all c G Consequently 


Jk+l — Jk+l9 — ‘^Jk9 Jk^gJkJ 

= 27fc - JfeCg*Cg7fc = 27fe - C4Cg7fc . 


The numerical computation yields fk=lk + Cfc, where Ck denotes the accumulated forward error. Let yk+i denote the new 
roundoff error, then in the k + I’th step of the iteration O we have 

fk+l — 2 fk Vk+l 

= lk+1 + (efc — Cg e/c) + (cfe — C*^ Cg jk) + Uk+i + C’dlefcH^) ■ 

The aforementioned estimates give 

\\i^ -C*k^9)^kh = IK-^- JkSg)ek \\2 < We-kh- 

Moreover the Janssen representation (see m p.l31]) gives 

C*ek ^9 Tfc = ^ X! Ti/p g)Mg/^ Ti/p 7fe. 

The last expression, when 7 ^ is replaced by 7 °, is the orthogonal projection Ify^efe of Cfc onto Vg := span(C/(p, 1//3,1/a)) 
and therefore 

\\ek-ClCg^°h<\\ekh 

and 

llefc - Cl Cg 7 fe ||2 < llefe - C^ Cg 7“||2 + \\Cl Cg ( 7 ° - lk)h < + 0{q^"). 

Since \\"f° — "fk \\2 = 0{q'^ ) by (|5]l, the new numerical error is 

||efc+i||2 < (1 + q^'^ekh + hk+ih + OiWekWl +9^'). 

In contrast to the operator version, ||efc ||2 enters linearly with coefficient si 1. Thus the numerical stability is plausible and is 
also confirmed by our numerical results in Example [3] 

(iii) Janssen’s alternative vector version: Janssen im Algorithm IV] proposes the approximations 

70 = A 5 , 

7fc+i = 2 7fcfc G No- (8) 

However, in their numerical tests of (l 8 ]l Janssen and Spndergaard ||9l observed some numerical instability. With notations as 
in (ii), the numerical error in step A: + 1 is 

Cfc+i = 7fc+i ~ 7fe+i 

= 2 Cfc — Cl Cek g — Cl C^y. g + yk+i + 0 (||efc|| 2 )- 

We show that the error may grow by at least a factor of 2 in each step. Note that Vg in (ii) is a proper subset of L^(K.), if 
a/3 < 1. Hence ||(/ —Hy^) ei ||2 ~ e, where e denotes the machine accuracy. The Janssen representation implies, that Ce^, g 
and Cl C^f, 5 are in V/. Therefore, 

(/ - Hvg) Cfe+i = 2 (/ - Hy^) Cfe + (7 - Vvlyk+i + C>(||efe||2) 

This shows, that the error components in can double in each step. In Example [3] we demonstrate that this numerical 
instability may indeed occur. 
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MONSTER canonical dual ol MONSTER 



Fig. 1: MONSTER function and its canonical dual for a = 20 and /3 = 1/50, computed by a routine of LTFAT IfTTI . 



Fig. 2: Error \\Cjo —CxfJ\\ of the three implementations of Algorithm |2] for 20 iteration steps to approximate the canonical dual 
window of MONSTER for a = 20 and /3 = 1/50. 


Example 3. As it was described in (iii) the implementation ® by Janssen can have some stability problems and should be 
applied carefully. We use the function MONSTER (see Figure [T]) in HI and a = 20, /3 = 1/50 for our numerical tests of all 
three implementations of the Schulz iteration. As we can see in Figure |2] the operator version is stable, the vector version of 
(ii) is also useful, while the error of the implementation in (iii) explodes. The same conclusions hold for the window function 
g{x) = e~'^^ with a = 20,/3 = 1/50, as is shown in Figured 



Fig. 3: \\C^° — of the three implementations of Algorithm |2] for 20 iteration steps to approximate the canonical dual 
window of the Gaussian g{x) = for a = 20 and /3 = 1/50. 
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II. Gabor frames of totally positive eunctions and exponential B-splines 


We now consider totally positive functions of finite type and exponential B-splines as window functions. The TP functions 
attracted much interest recently, as they provide new examples of window functions for which the necessary density condition 
a/3 < 1 of the lattice parameters is also sufficient IfTSl . For detailed information on total positivity of functions and matrices 
see OH, for a detailed introduction to exponential B-splines see ll20i and ED. 

Definition 4. Il22l . Il23l An integrable function g : K —?> R is called totally positive (TP), if its Fourier transform factors as 




/ dt 

J —OO 


°° ^2ni6t,uj 




1 + 2TTiS,^ui ’ 


(9) 


where C, rj, 6, are real parameters with 


OO 

C > 0, r] > 0, 0 < rj + ^ Sl < oo . 

i/=i 


Note that in Schoenberg’s terminology there also exist TP functions that are not integrable. Therefore the given definition 
does not include the most general case of TP functions. Here we only consider TP functions of finite type N G N, with 
ry = 0, A = (i5i,..., 6n) and 

N 

+ 2m5vOj)~^. ( 10 ) 

v=\ 


For N = 1 and A = {6) we have the one-sided exponential 

55 ( 2 ;) = '"^X(o.oo)(fe) for a::eR\{0}, 

with support [0, 00 ) for positive S and (—oo,0] for negative S, and for A > 1 and A = {Si,... ,6 n), Pa is the A-fold 
convolution 

9A = gsi * ■ ■ ■ * gsN- 

Especially, if 6,^ = S 0 for all 1 < i/ < N, then 

U|N -1 _ 

The functions pA are nonnegative, have infinite support, and exponential decay, precisely, if t = (max{ 27r |5jy| | 1 / = 
1,..., A})“^ and e > 0, then there is a constant c = c(e) such that 

0<PA(a;) (11) 


Moreover, if 6n and we specify p 5 ( 0 ) = 1/|25|, the recurrence relation 


9Si,...,Sn = 




- s 


N 


holds. For later use, we denote by m (resp. n) the number of positive (resp. negative) parameters 6,^. 

The main result in lITSi shows that the Gabor system Q{g, a, /3) of a TP function of finite type m + n > 2 constitutes a Gabor 
frame if and only if a/3 < 1. The proof provides an algorithm for the computation of a dual window 7 with compact support. 
This method was adapted in ll24l . Il25l in order to supply infinitely many dual windows 71 ,. The computation of 7 ^( 2 : + ja), 
with j G Z, is performed by computing a single row of a left-inverse of the biinfinite pre-Gramian matrix 


P,(x) + 


( 12 ) 


The structure of the left-inverse of Pg{x) heavily depends on the property that g is totally positive. We include the algorithm 
for the reader’s convenience. 


Algorithm 5. II24I Input parameters are the parameter vector A = (Si, S 2 , ■ ■ ■, Sjv) of the window g, the lattice parameters 
a,/3 > 0 with a/3 < 1, a parameter L gNq controlling the support size of the dual window 7 l, and a point x G [0,a). 

Output parameters are integers ii{L), 12 (T) and the vector of values 7 ^( 0 ; + cvj), ii{L) < j < i 2 {L), in the support of 
such that 7 l : M —>■ R. defines a dual window of g. 


1) Set r := 


1 —a/3 


J, fci = —(r + f)m and k 2 = (r + l)n. 
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2) Set ki{L) := ki — L and k 2 {L) := k 2 + L, 

ii{L) := 

i2{L) := 


ki — L + m — 1 

al3 

k2 L — Tl 1 


a/3 


X 

a 

X 

a 


+ 1 ; 

- 1 . 


3) Set Pl{x) := {Pj,k)ii{L)<j<i 2 {L), ki{L)<k<k 2 {L), where 


Pj,k = g[x + aj - 

4) Compute the pseudoinverse 

Pl{x)‘^ = (9fcj)fci(L)<fc<fc2(L), ii{L)<j<i2{L) 

of Pl{x). 

5) Take the row with index k = 0 of PL{x)f Its coefficients define the values of the dual window 7 l at the points 
{a; + aj \ ii{L) < j < i 2 {L)}, i.e. 


')l{x + aj) := 


/3go.i , if h{L) < 3 <i2{L), 

0 , if j < ii{L) or j > i 2 {L). 


(13) 


In particular, the support of the dual window 7 ^ is contained in the interval [azi(L),a(i 2 (i) + 1)] of length of the order 
(3~^{k2 — ki + 2L). Thus the parameter L labels the size of the support. 

A related class of window functions is the class of exponential B-splines. These functions are positive and have compact 
support, a property which is desirable in some applications. 

Definition 6. For A = (Ai,..., Xn) € the exponential B-spline (EB-spline) with knots 0,1,..., is given by its 
Fourier transform 


N 


bk { u 2 )=n 


I/=l 


pAix—27riu; _ ^ 

Xu — 2'Kiui 


(14) 


In Qa, it was shown that the Gabor system G{Ba, a, /3) of every EB-spline constitutes a frame for a = 1, /3 < 1 (and also 
some other lattice parameters). Similar to the case of TP functions, the following algorithm provides dual windows 7 ^ of the 
Gabor frame C/(i?A, a,/3). 

Algorithm 7. lHa Input parameters are the parameter vector A S of the window Ba, the lattice parameters a, /3 > 0, 
a parameter L gNq controlling the support size of the dual window jl, ond a point x S 

Output parameters are integers ii{L), i 2 {L) and the vector of values ^l{x + aj), ii{L) < j < i 2 {L), in the support of 
such that 7 a : M —>■ K. defines a dual window of Ba- 

1) Set k 2 {L) = j + I + L and ki{L) = —k 2 {L). 

2) Set ii{L) < ki{L) and i 2 {L) > k 2 {L), such that 

Ba [x + (^l(T) - l)a - ^ q, 

SA(x + ii(L)a-^i^) = 0 , 

Ba [x + ML) + l)a - ^ q, 

BA(x + i2{L)a-^^^^) =0. 

3 ) Set Pl{x) := iPj,k)ii{L)<j<i2{L), ki(L)<k<k2(L)> where 

k 


Pj,k = Ba^xX- aj - - j . 

4) Compute the pseudoinverse 

Pl{x)^ = (gfc,j)fci(L)<fc<fc2(i), ii{L)<j<i2{L) 

of Pl{x). 

5) Take the row with index k = 0 of PL{x)f Its coefficients define the values of the dual window 7 l at the points 
{a; + aj \ ii{L) < j < * 2 (A)}, he. 


7 L(a; + aj) := 


Pqo,j , if h[L) < j <i2{L), 

0 , if j < ii{L) or j > i 2 {L). 


(15) 
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We would like to emphasize the following point: (i) These algorithms determine the precise values of a dual window 7 
with compact support in L^(R) and not just a discrete approximation in a finite-dimensional vector space, as is done in most 
existing algorithms in Col, im. 

(ii) Although the problem of finding a dual window is by nature infinite-dimensional in L^(K.), the computation of on 
a grid requires only the pseudo-inversion of M finite-dimensional matrices. 


III. Approximation of the canonical dual of TP functions and EB-splines 

We fix the parameters a, (3 > 0 of the Gabor frame and let g be either a TP function of finite type with parameter set A 
or an EB-spline with parameter set A. The pre-Gramian matrix (fT2l) plays a central role in the characterization of the frame 
bounds A,B of the Gabor system Q{g,a,l3), and in finding dual windows 7 . 

In this section we address the question how the sequence of dual windows 7 ^ of Algorithms |5] and [T] are related to 
the canonical dual window 7 °. Our main goal is to prove that the compactly supported dual windows computed by the 
Algorithms |5] and I 2 ] from ll24l and ifThl approximate the canonical dual window 7 ° at a rate 

hi -7°ll2 < ce-P^, 

where L determines the support of 7 ^ and p > 0. 

As a first step we show that the norm of 7 ^ remains bounded as L tends to 00 . 

Theorem 8. Let g be a TP function of finite type as defined by ( 1701) . Then there exist constants A, B > 0 independent of 
L € No such that 

>l||c ||2 < ||7^L(a;)c||2 < 5 ||c||2 (16) 

for all c € ^ G [0, a), where Pl{x) is the finite section of the corresponding biinfinite pre-Gramian matrix 

Pg{x) as described in Algorithm]^ Consequently, ||7l||2 A y/a A~^, for all L G N. 

The precise proof will be given in the appendix. The proof idea goes as follows: The upper bound B is easily obtained 
from Schur’s test based on the exponential decay of g in (fTTl l. Eor the lower bound A, we choose 2L -\-1 submatrices Ri of 
Pl{x), with 7:2 — fci + 1 columns each, and build a left-inverse Ql{x) of Pl{x) by the selection of specific rows of 7?J. The 
result in ifTSl Theorem 9] shows that all matrices are bounded uniformly in x. Their upper bound C will provide the upper 
bound A~^ = C for all matrices Ql(x) uniformly in x and L. Hence A is a suitable constant for the lower bound 

of Pl{x). 

To prove the rate of approximation of the dual windows 7 ^, we recall that the canonical dual window 7 ° = Sf^g can be 
expressed by the Moore-Penrose pseudoinverse of Pg{x), namely 

7 °(a: -f aj) = P (Pg(x)’l')oj (17) 

for all j G Z and all x G [0, a). This is a direct consequence of the Wexler-Raz criterion for the dual windows of p || 8 ] 
Theorem 7.3.1], and the minimal L^-norm of the canonical dual among all duals of g EEj . We will therefore show that 
the zeroth row of PLix)^ approximates the zeroth row of Pg{xY at an exponential rate. Eor this we need a new result on 
non-symmetric finite sections of biinfinite matrices. The following theorem is not covered by the results in Il27l and may be 
of independent interest. To fix the notation, for n = (ni,n2) G N^ and b G ('^(Z), we let Vn with 

Pr^b^f. ., 0 , ,..., bYi2—i, bji2 5 O 7 • ■ •) 

be the orthogonal projection onto the ni n 2 + 1 -dimensional subspace 7^„f^(Z) = C"i+«- 2 -i-i Pqj- ^ biinfinite matrix 
U = {uj^k)j,k& and r,n G N^, Vr UPn is a non-symmetric finite section of U. We will write {Pn UPn)~^ for the inverse 
of the symmetric finite section Pn U Pn on the finite-dimensional subspace Pn P' (Z) (with the understanding that it cannot 
be invertible on f^(Z)). 

Theorem 9. Let {x{P)) ® strictly increasing sequence of integers and U = {uj^k)j,k& be a biinfinite matrix such that 

(a) U*U is invertible, and (b) there exist constants c, a > 0 such that 

\u 3 ,k\ < for all j,kGZ. (18) 

Let / C and assume that for every n G I a finite section Un '■= Pr{n) U Pn A given such that 

^4 ||c ||2 < ||(7„ CII 2 for all c G ilj) (19) 

for some constant A > 0 independent of n. Then there are constants c, d > 0, such that for all n G I 


\\U{U* U)-^ eo - Un {U*n Un)-^ eo||2 < £6-“”“ 


( 20 ) 
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where uq := mm{ni,n 2 ,ri(n),r 2 (n)}. 

The proof is deferred to the appendix. 

Note that the row vector [U{U*U)~^)eo)* = Cq {U*U)~^U* = is precisely the zeroth row of the Moore-Penrose 

pseudoinverse of U as it arises in the computation (fTTI l of the dual windows 7 l and 7 °. We also note that the decay condition (fTST l 
models the decay of the entries off a “ridge” x(^) rather than off-diagonal decay (in which case x(^) = ^)- The above 
definition reflects exactly the behavior of the pre-Gramian matrix Pg{x) of a window with exponential decay, since \Pg{x)jk\ = 
\g{x + aj — k//3)\ < in which case x(^) = [k/{aP)\- Decay conditions of this type occur in wavelet 

theory ll28l and in the theory of Fourier integral operators ll29l . 

We can now prove the main result of our paper, namely that the numerically computable dual windows 7 ^ converge 
exponentially fast to the canonical dual window 7 °. 

Theorem 10. The dual windows 7 ^, L G N, in Algorithm\^approximate the canonical dual window 7 ° of g at an exponential 
rate 

\hL-l°\\2<ce~P^. ( 21 ) 

Proof: We set n{L) = (|fci|-|-L, fc 2 -l-T) and r(L) = {\ii{L)\,i 2 {L)) as in Algorithm |5] The matrices (cc) in Algorithm|5] 
are exactly the non-symmetric finite sections Pl{x) = Pr(L)Pg{x)Pii(L) of the pre-Gramian Pg{x). Then Theoremimplies 
that the finite sections Pl{x) of the pre-Gramian Pg{x) are left-invertible (on the appropriate finite-dimensional subspaces) 
with constants independent of L. Therefore the decay conditions and the uniform bounds in the assumptions of Theorem 0 are 
fulfilled. Thus for fixed x we obtain that 


Finally the approximation of the dual windows 7^ and 7° follows from 

/ OO 

\il{x) - "i°[x)\^ dx 

-00 

/*a 

= / \il{x + aj) - 7°(x -f aj)f dx 
Jo 


= 13'^ [ WeoPUxy - e^Pg{x)'<\\ldx 
Jo 


Since no(L) := min{ni(P), n 2 (L),ri(L), r 2 (L)} = min{|fci| + L,k 2 + L, \ii{L)\,i 2 {L)} = L + C for some integer constant 
depending on the window only, the rate of approximation in (I 2 TI 1 follows. ■ 

Remark 11. (i) Theorem |9] is not contained in the results on the non-symmetric finite section method in Il27ll . The selection of 
rows by r(n) in our assumption (fT9] l meets only the condition in lIZTl Lemma 5.2], which reads as 

sup||(P* < OO 

nGl 

in our notation. However, the condition in lIZTl Lemma 5.1] is not matched and can only be satisfied by considerably increasing 
the number of rows of The approximations of the canonical dual 7 ° based on the non-symmetric finite section method in 
ll27l do not provide dual windows, in contrast to our approximations jl. 

(ii) With an analogous proof, we obtain that the dual windows 7 ^ in Algorithm [T] approximate the canonical dual of the 
EB-spline Pa at an exponential rate. 


IV. Discretization and implementation 

For the numerical use of TP functions and EB-splines in signal analysis it is often necessary to discretize these windows. 
We define the sampling operator Ss for a given sampling rate 6 > Oby 

Ssg := Vs {g{6k))^^^ 

and the periodization operator Pr- with period K > Oby 

PKg{x) ■-'^g{x + Kk), xG[0,K). 


fcGZ 


Eor discrete signals c G V (Z) we let 


pkc-.= { y : 


Cj-\-Kk I € C 

J j=0,...,K-l 


K 
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Furthermore we consider the dilation operator 

Dhg{x) := Vhg{hx), h > 0, 

which preserves the L^-norm. The exponential decay (HB of TP functions or the compact support of EB-splines imply that Ss g 
is in C £^(Z) and Pxg G L^([0,iT)). The combination of both operators yields finite discrete signals PxSsg G C^. 

For a, b, M, N,K G N and Mb = Na = K, the discrete Gabor system is defined as 

g{PK Ssg,a,^) = Pk Sg g{- - ka) 

The following result is well known and holds for an extremely general class of window functions. All TP functions of finite 
type and all EB splines together with the duals of Algorithms |5] and |7] satisfy this mild condition. 

Proposition 12. Let a, /3 > 0 and a/3 = ^ = -^ o.nd Mb = Na = K with a, b, M, N,K G'H. Let 5,76 such that 

{g(x + ja))jgz and {^{x + ja))j^z are absolutely summable for all x G [0,a), and 

^y{x + ja)g{x + ja - kjf) = /3i5o,fc 

jGZ 

for all x G [0, a) and fc G Z. Then g{PK S^/a d: a, is a Gabor frame for and ^{Pk S^/a 7i a, ■^) is a dual Gabor 
frame. 

This statement is proved under slightly different assumptions in j^. The assumptions in Proposition [12] lead directly to 
the verification of the Wexler-Raz criterion in IMl Theorem A.3] for dual Gabor frames of C^. The details of the proof are 
omitted here. 

Remark 13. Since SgDh = Sgh, it is helpful to dilate the function g by the sampling rate. Subsequently we can work with a 
sampling rate (5 = 1 and consider PxSig =: PrScj of some scaled TP function or EB-spline g := D^/i^g. In many practical 
situations a/a is proportional to s/K. Thereby the time-frequency localization of the window is independent of K. 

In the remaining part of this section we describe implementations of discretized TP functions and EB-splines as well as 
the duals from Algorithms |5] and [Tj in Section In] For this purpose, we use some knowledge about the Zak transform of these 
functions. For a parameter a > 0 and a function / G L^(R) with absolutely summable {f{x -I- aj))^^^ for a: G M, the Zak 
transform is defined by 

Z„/(x,a;) :=^/(x + aj)e-2™. 

jGZ 

The Zak transform is a-quasiperiodic in x and 1/a-periodic in w m- For a given periodization parameter AT G N we obtain 
the discrete version of a scaled TP function or EB-spline g by 

PKS~g{k) = ZKgik,0), k = 0,...,K-l. (22) 


k = 0,..., N — 1 and 
I = 0,...,M -1 


A. EB-splines 

Since EB-splines have compact support, their Zak transform is a finite sum and only requires finitely many point evaluations 
of these functions. 

Case 1: The EB-spline \ with a single weight A G K of multiplicity m G N can be factorized into an exponential and 
the cardinal polynomial B-spline Nm of order m 

= (X[o.i) * • • • * X[o,i)) = 

Hence it can be evaluated by the well-known algorithm by Cox and deBoor (see JIB), which is part of the standard signal 
processing toolboxes. 

Case 2: For EB-splines with pairwise distinct weights Ai < • • • < Am and m > 2 Christensen and Massopust |[33ll give the 
closed form 

m 

Bxi,...,Xmi^ -f fc — 1) = ^ x G [0,1), 1 < A: < m, 

i=i 


with coefficients 


n ~ ^r) , 


fc = 1, 


(k) 

= < 


= 1 , 

E 

31W 


•^71 + • ■ • + 7 . 


r=l, 


fc = 2,..., m. 
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Case 3: For EB-splines with several distinct weights with multiplicities, we use the following four-term recurrence relation, 
which was stated in ll34l and ll35l . 

Theorem 14. Il34l . Il35l Let Ai,..., Aat S K, Xi ^ Xn and 


the associated EB-spline. Then the following recursion holds 

/ \ ^Ai,...,AAr-i (^) -^A2,...,Aiv (^) , ^ ^^A2 ,...,Aa;- 

Ai — An 


1 ) - 
Ai — Xn 


1 ) 


By iteration of this recurrence it is possible to reduce any given EB-spline into either several lower order EB-splines 
with pairwise distinct weights of multiplicity one or only one weight of higher multiplicity. These can be treated as in the 
aforementioned cases. 


B. TP functions 

In the case of TP functions, we present two different implementations of the computation of PxSg in (l22ll . Eor both we 
use that TP functions are invariant under dilation. 

Lemma 15. Let g be a TP function of finite type N gN with weights and h > 0 a scaling parameter. Then 

gh = VhDhg 

is the TP function of finite type in (doll with weights 

The first implementation uses the identity in ll24l Remark 2] 

7jKgh{x,w) = ..., I Vx e [0,iT), 

where the right-hand side is the divided difference of with 

( N \ —hxy 

1 fKihy+2.^.) 

in the knots ^,..., . 

The second implementation uses the connection of TP functions to EB-splines. In ifThl Theorem 3.4] it is shown that the 
Zak transform of a TP function can be expressed in terms of the Zak transform of an associated EB-spline. 

Theorem 16. |1I§| Let g be a TP function of finite type with weights Si,... ,6^ G R- With A^ := — i/ = 1,..., N, we 
have 

KZKgh{x,0) = l[ ^^^Zii3A„....A„(f,0), xG[0,1). 

6 1 

iy—1 

Consequently PrSp can be computed by the Zak transform of the corresponding EB-spline as described in IIV-AI 
C. Dual windows 

Eor TP functions g and lattice parameters a, (3 > 0 with aP < 1, Algorithm |5] in Section HJ allows us to compute samples 
Sa/NlL, G N, of a dual window jl. Likewise, Algorithm|7]provides the sampled dual windows of EB-splines B\. Therefore, 
for a given periodization parameter AT G N and the time-shift parameter a G N, we compute 

PKSa/alL{k) = E,{Sa/al L)[k + j K), fc = 0, . . . , AT - 1 . 

jGZ 

This sum is finite because of the compact support of 7 l. 

Example 17. We use Algorithm |5] for the computation of the dual Gabor window of the asymmetric TP function g with 
parameters 5 = [—1,1,1/3,1/5] and lattice parameters a = 2/3, P = 1. The discretization parameters K = 900, a = 20, 
b = 1/30, are chosen according to the standard dilation aja = s/K in Remark fTH Eigure |4] shows the discrete TP function 
PK^a/ag and its dual window PRSa/alL for L = 20. The difference \\PKSa/a{lL — 7°)||2 to the discrete canonical dual is 
7 • 10“® measured in the .f 2 -norm of 


Sometimes the frame bounds in the discrete case may be better than in the continuous case. Therefore the discrete TP 
functions or EB-splines may even provide Riesz bases at the critical density, as is explained in the following recent result. 
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TP function g 




Fig. 4; Discrete TP function and its dual 720 for a = 20 and M = 30. 


Proposition 18. Il24l . ll3^ Let g be a continuous TP function as in including the infinite type with rj = 0. Assume 
a = M G N and let = 1/M and Ff £ N such that K/M € N. If K/M is odd, then Q{PKSg,a, fi) is a basis ofC^. 

In addition, assume that g is even, which means that {5^ | > 0} = {—5^ \ by < 0}. If M is odd, then Q{PKSg,a, fi) is 

a basis of C^. 

Remark 19. In the critical case afi = 1 the computation of the dual window 7 cannot be performed by the aforementioned 
discretization procedure, as the Gabor system Q{g,a,l3) is not a frame in L^(R). In this case, the usual method using the 
discrete Fourier transform and the discrete Zak transform of PxSg should be applied for the computation of the discrete dual 
window EH. 
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Appendix 

In the appendix we provide the technical details of the proofs of the main theorems. The proof of Theorem 0 depends 
heavily on ifTSl Theorem 9]. 

Proof of Theorem^ The pre-Gramian is bounded on as a consequence of the exponential decay of g and Schur’s 

test. Therefore the finite sections Pl{x) = Pr(n) Pg{x)Pn are uniformly bounded. 

To prove the existence of a lower bound A, we construct a left inverse of Pl{x) by adapting the proof of lITSl Theorem 9]. 
By step 2 of Algorithm |5]FL(a;) has columns indexed hy k, ki — L < k < k 2 + L and rows indexed by j, ii{L) < j < i 2 {L), 
thus every left-inverse Ql{x) has columns indexed by ii{L) < j < i 2 {L) and rows indexed hy ki — L < k < k 2 + L. We 
will construct these rows one by one in three steps. 


5tep 1. For every index —L < i < L, we choose the following submatrix of Pl{x). Given x and f £ Z, we choose 
ye £ [ 0 , a) and je £ Z, such that 

ye = x + jea - 

Then we define the matrix Re = Po{ye) by using the same definitions as in step 3 of Algorithm |5] with L = 0; i.e., we let 


*i(0) 

4(0) 


fci 4- w - 1 _ ^ 
a/3 a 

k2-n+l _ ^ 
a/3 a 


and 


Re (0)<4<i2(^) J k\<,k<k2 ’ 


ye + ja 


k . . I Pk 

- = X + [j + je)a - — 


with re-j^k = 9 [ye A- ja - | j. Note that 
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Hence, the columns of Ri, as indexed by fci < fc < ^ 2 , correspond to sections of columns of Pl{x) indexed by ^ + fci < fc < 
£ + k 2 - More precisely, R-l contains sections of the brst k 2 — ki + 1 columns of Pl{x)', with increasing £ the selection of 
columns is shifted to the right; and Rl contains sections of the last k 2 — ki + 1 columns of Pl{x). For the rows of Ri, we 
observe that 

+ l = i{{0)+ji 

and i 2 {L) > 1^(0) + je- To summarize, each Rg is a ( 12 ( 0 ) — i{{0) + 1) x {k 2 — fci + l)-submatrix of Pl{x) with a dimension 
independent of \£\ < L. 


ii{L) < 


ki + m — 1 X I 

a/3 a a/3 


5tep 2. The arguments in IfTSl Theorem 8] relate left-inverses of Ri = Po{ye) to left-inverses of the biinhnite matrix Pg{yi). 
It is shown that 

• Ri has full column rank, 

• there exists a uniform bound C > 0, which does not depend on yi G [0,a), and left-inverses F^ of R^ such that 

IIF^II < C for all - L < i < L, 


• the rows with index fci < /c < 0 of F^ are orthogonal to all columns with index k' > k2 — n of the biinhnite matrix 
Pg{yi), and hence orthogonal to all columns k' > l + k2 — n of Pl{x). Likewise, rows 0 < fc < A:2 of F^ are orthogonal 
to all columns with index k' < £ + ki + m of Pl{x). 

5tep 3. With these properties, we obtain the left-inverse Ql{x) of Pl{x) by dehning the rows of Ql{x) as follows: 

• We start with r_L and take its rows fci < fc < 0 as the hrst rows of Ql{x), extended by zeros such that 


{QL{x))k-L,j+j-L ^ 



*1 {0)<j<i2 ( 0 )) 

Otherwise. 


• For —L + 1 < £ < L — 1 we take the row with index 0 of F^ and extend this row by zeroes. 


{QL{x))ej+ji 


0, otherwise. 


• We end with F l and take its rows 0 < fc < fc 2 as the last rows of Ql{x), extended by zeros such that 


{QL{x))k+L,j+jL 


(rL)fcj, i{'(0) < j < i^(0), 
0, otherwise. 


It is clear that all entries of Ql{x) are bounded by the constant C in step 2. Moreover, every row and column of Ql{x) has 
at most 

fc2 — fci + 1 


max (info) — iUO) + 1 ) < 
-L<l<L 


aP 


nonzero entries. Therefore, by Schur’s test, we obtain that 

\\QLix)\\<C 


k2 — ki + 1 


Since is a left inverse of Pl, we obtain 

||c ||2 = \\Ql{x)Pl{x)c \\2 < c 


a/3 

fc2 — fci + 1 
a/3 


ll^’L(a;)c||2, 


and we may choose A = {C ^ as a lower bound in (fThl l. This completes the proof of (fThl) . We conclude that 

||eJPL(x)t ||2 < =24-1. Therefore 


/ OO 

\lL{x)\^dx 

■00 

= P^ [ '^\il{x + dx 
•10 jGZ 

= /3^y \\eQ Pl{x)''^ Wl dx < a A~'^. 


In the following we will need a well known result of Jaffard Il3^ and Baskakov ll^ . 
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Proposition 20. Let A = 

C,a>0, such that 

If A is invertible, then A~^ 


{Ajk)j^k& be a biinfinite matrix with exponential off-diagonal decay, i.e., there exist constants 


yj,k€Z. (23) 

also possesses exponential ojf-diagonal decay and there exist C" > 0 and 0 < a < a, such that 
|(A-i),fe| < Vj,fcGZ. 


Corollary 21. Let A be positive and invertible on £^(Z) with exponential ojf-diagonal decay (I231 l and let VnAVn be a 
sequence of finite sections. Then every matrix VnAVn is invertible on and there exist C > 0 and d,0 < d < a, 

such that 

I {VnAVn)-k I < for -n,<j,k<n 2 . (24) 

Proof: The proof follows BSl and Eli. Let A be the matrix obtained by stacking the finite sections VnAVn along 
the diagonal. Then A possesses exponential off-diagonal decay (|2^ . Since A is invertible and positive, its spectrum a{A) 
is contained in an interval [^,5] for some A,B > Q, consequently the spectrum of the restriction of VnAVn on P„£^(Z) 
is also contained in [A, B] and every finite section is invertible on VnP'i'^)- Therefore the stacked matrix A is invertible 
on 07^„('^(Z) ~ By Proposition l20l A~^ possesses exponential off-diagonal decay. Since A~^ consists of the blocks 

{VnAVn)~^, they satisfy \{VnAVn)~k \ < for -ni < j,k < n 2 . ■ 

We remark that clearly every finite matrix possesses exponential off-diagonal decay. The point is that the constants may be 
chosen independently of the size of the finite section. 

Proof of Theorem |9} Recall that, for n = (ni, 712 ) G N^, Vn b = (..., 0, ..., 0,.. .)^ is the 

orthogonal projection onto 7^„£^(Z) = C”i+” 2 +i We write Un = VrUVn for a non-symmetric finite section of U. In the 
assumption of Theorem depends on n, but we will omit this dependence in the notation. All operator norms are in f^(Z). 
Let b := U {U* U)~^ eo and dn ■= Un {Ui^ Un)~^ cq. We decompose the norm into three parts 


\\b-dnh<\\{U-UVn){U*U)- 


eo 


-f 


UVn (iU* [/)-! - {Vn U* UVn)-^) Bq 


P\\UVn{{VnU*UVn)-^ -Un{U*M~^)^A 


2 ■ 


(25) 


1 . Note that Lf* U is invertible and positive. Since U fulfills the decay property (fT^ . it is easy to see that the symmetric matrix 
U* U decays exponentially off the diagonal, i.e., for some constants C,a > 0 

|([/*L)j-fc| < j,k€Z. 

By Proposition l20l the inverse matrix inherits the exponential decay, and thus there exist constants C" > 0 and a,0 < a < a, 
such that \ {U* U)J^\ < for all j, fc G Z. Hence the entries of the vector v := {U* U)~^ eo also decay exponentially 

as 

kil < Cl j G Z. (26) 

Since 

— np —1 00 

j=-oo j=n2 + l \j\>no 

the decay property (l26l l implies that 

\\{U-UVn){U*U)-^eo\\^ < ||L|| ||(/-iP,^)p||2 <C2e-“"L 


2. Since U* U is invertible and positive with spectrum a{U*U) C [A,i3] for A,B>0, the spectrum of the finite sections 
Vn U* U Vn on is also contained in [A, B], As in the finite section method in HOl and with v = {U* U)~^ eo, we 

obtain 

||((C/* U)-^ - {Vn U* UVn)-^) eoL = Wn U* U VnY^ Vn U* U {Vn " I) v\\^ < \\{U* U)-^\\ \\U* U\\ \\{Vn " I) t^|l2 

and consequently 

\\UVn {{U* U)-^ - {Vn U* UVn)Y < ^3 6 "“"“ 
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3. To treat the third term in (l25T l. we need a geometric interpretation of the rows of the Moore-Penrose pseudoinverse. Let 

br. ■■= {UVr.){Vr.U*UVr,)-^eo 

and dn as above. Then is the transpose of the zeroth row of the Moore-Penrose pseudoinverse of U Vn- By Corollary 1211 
i'Pn U* U'Pn)~^ satisfies (l24li independently of n. Therefore the same argument as in the first part implies the decay property 

|(&r^),| <C4e-“l^l. (27) 

It is essential that the constants are independent of n. Since the Moore-Penrose pseudoinverse of U Vn is also a left-inverse, 
we have 

bn -L Vn := spanjufc | -ni < k < n 2 , k^ 0}, 
where Uk, fc G Z, are the columns of the matrix U. Likewise for Un = VrUVn we have 

dn-LVrVn ■ 


Using this orthogonality, we rewrite the vectors bn and dn as follows. Let Hw denote the orthogonal projection onto some 
subspace W. Now set 

&„:=(/- ny„) uo and dn := (I - ^Vr.v„)VrUo. 

Since bn S lm{UVn) and dn G Im((7n), we obtain 


bn — 


and dn = 


The normalization in (l28l l is obtained from 

{bn 1 tTo) — 1 and (bnt t^o) — {bn j bn T tro) — II|| 2 ? 

{dn 5 V rtio) — 1 and {dn , V r Uq) — {dn 5 dn r tto) — || dn || 2 • 

As the third term in (l25T l equals ||6„ — d„|| 2 , we first consider ||6„ — d„|| 2 - For this purpose we write 

bn dn — (T V (tlQ Fly^ tlo) Vr Uq Vr Fly^ Uq 

= {I — Vr) bn + bl-p^Vr. Vr [uq — Uq) 

= (/ — Vr) bn + bl-p^Vr. 'Pr bn- 


(28) 


(29) 


By the assumption (fT9] l. the truncated columns Vr Uk, with —ni < k < n 2 and k 0, form a Riesz basis for VrVn with 
lower Riesz bound A. Furthermore it holds that 

2 


'Pr br 


<A ^ ^2 \{n'p^VnVrbn,VrUk)\^ 

— ni</c<n2, k^O 

= A-2 ^ \{Vrbn,Uk)\'^ 

— ni</c<n2, k^O 

= A-2 ^ \{{Vr-I)bn,Uk)\^ 


— ni</c<n2, k^O 


< A-^B 


{Vr - I) k 


where B = ||t7|p denotes the Bessel bound of all columns Uk, fc G Z. Taking the f^-norm in (l29T l and substituting the above 
estimate, we obtain 

Wbn-dnh < C5\\{I -Vr)bn\\2- 

We now return to \\bn — dn\\ 2 - B is an easy exercise that for nonzero vectors y,w we have 


y 

W 

/ 3 || 2 /- 

W\\2 

Ml 

Ikll2 

2 ~ min{l|y|l2 

,|Hl2} 


Using once more that bn G lm{UVn) and dn G Im([4i), the assumption (fT9l l implies that 

0 < A < ll^nlb, Wdnh < lko|l 2 ■ 

By the decay property (l27l i we obtain 

||&^-dr.||2 < ^||L-dr.||2 < ^ || (/" Pr ) L || 2 < Cg 6"“. 

To finish the proof of Theorem |3 we add the three contributions in (|25]) and obtain the convergence rate (l20l) with a constant 

C = C 2 + C 3 -I- C 6 . ■ 



















16 


References 

[1] M. Dorfler, “Time-frequency analysis for music signals: A mathematical approach,’’ J. New Music Res., vol. 30, pp. 3-12, 2001. 

[2] R C. Russel, J. Cosgrave, D. Tomtsis, A. Vourdas, L. Stergioulas, and G. R. Jones, “Extraction of information from acoustic vibration signals using 
Gabor transform type devices,” Meas. Sci. TechnoL, vol. 9, pp. 1282-1290, 1998. 

[3] F. Hlawatsch and G. Matz, Wireless Communications over Rapidly Time-Varying Cannels. Academic Press, 2011. 

[4] S. Sriram, N. Vijayakumar, P. A. Kumar, A. S. Shetty, V. P Prasshant, and A. Narayanankutty, “Spectrally efficient multi-carrier modulation using Gabor 
transform,” Wireless Engineering and Technology, vol. 4, pp. 112-116, 2013. 

[5] S. Blanco, C. E. D’Attellisand S. I. Isaacson, O. Rosso, and R. O. Sirne, “Time-frequency analysis of electroencephalogram series, II: Gabor and wavelet 
transforms,” Physical Rev. E., vol. 54, pp. 6661-6672, 1996. 

[6] L. Chen, E. Zhao, D. Wang, Z. Han, S. Zhang, and C. Xu, “Feature extraction of EEG signals from epilepsy patients based on Gabor transform and 
EMD decomposition,” in Sixth International Conference on Natural Computation, 2010. 

[7] O. Christensen, An Introduction to Frames and Riesz Bases. Birkhauser, 2003. 

[8] K. Grochenig, Foundations of Time-Frequency Analysis. Birkhauser, 2001. 

[9] A. J. E. M. Janssen and P L. Spndergaard, “Iterative algorithms to approximate canonical Gabor windows: computational aspects,” J. Fourier Anal. 
AppL, vol. 13, pp. 211-241, 2007. 

[10] T. Strohmer, “Numerical algorithms for discrete Gabor expansions,” in Gabor Analysis and Algorithms, H. G. Feichtinger and T. Strohmer, Eds. 
Birkhauser, 1998, pp. 267-294. 

[11] Z. Prusa, P L. Spndergaard, N. Holighaus, C. Wiesmeyr, and P Balazs, “The large time-frequency analysis toolbox 2.0,” in Sound, Music, and Motion, 
ser. Lecture Notes in Computer Science, M. Aramaki, O. Derrien, R. Kronland Martinet, and S. Ystad, Eds. Springer International Publishing, 2014, 
pp. 419^42. 

[12] G. Schulz, “Iterative Berechnung der reziproken Matrix,” Z. Angew. Math. Mech., vol. 13, pp. 57-59, 1933. 

[13] A. J. E. M. Janssen, “Analysis of some fast algorithms to compute canonical windows for Gabor frames,” 2002, unpublished. 

[14] -, “Some iterative algorithms to compute canonical windows for Gabor frames,” Proceedings of IMS Workshop on Time-Frequenzy Analysis and 

Applications, 2004. 

[15] K. Grochenig and J. Stockier, “Gabor frames of totally positive functions,” Duke Math. J., vol. 162, pp. 1003-1031, 2013. 

[16] T. Kloos and J. Stockier, “Zak transforms and Gabor frames of totally positive functions and exponential B-splines,” J. Approx. Theory, vol. 184, pp. 
209-237, 2014. 

[17] K. Grochenig, “Acceleration of the frame algorithm,” IEEE Trans. Sign. Proc., vol. 41, pp. 3331-3340, 1993. 

[18] H. Hotelling, “Some new methods in matrix calculation,” Ann. Math. Statist., vol. 14, pp. 1-34, 1943. 

[19] S. Karlin, Total Positivity. Stanford University Press Inc., 1968. 

[20] L. L. Schumaker, Spline Functions: Basic Theory. John Wiley and Sons, 1981. 

[21] R. A. DeVore and G. G. Lorentz, Constructive Approximation. Springer, 1993. 

[22] I. J. Schoenberg, “On totally positive functions, Laplace integrals and entire functions of the Laguerre-Polya-Schur type,” Proc. Natl. Acad. Sci. USA, 
vol. 33, pp. 11-17, 1947. 

[23] -, “On Polya frequency functions. 1. The totally positive functions and their Laplace transforms,” J. Anal. Math., vol. 1, pp. 331-374, 1951. 

[24] S. Bannert, K. Grochenig, and J. Stockier, “Discretized Gabor frames of totally positive functions,” IEEE Trans. Inform. Theory, vol. 60, pp. 159-169, 
2014. 

[25] T. Kloos, “Gabor Frames total-positiver Funktionen endlicher Ordnung,” Diploma thesis, TU Dortmund, 2012. 

[26] A. J. E. M. Janssen, “Duality and biorthogonality for Weyl-Heisenberg frames,” J. Fourier Anal. AppL, vol. 1, pp. 403^36, 1995. 

[27] K. Grochenig, Z. Rzeszotnik, and T. Strohmer, “Convergence analysis of the finite section method and Banach algebras of matrices,” Integr. Equ. Open 
Theory, vol. 67, pp. 183-202, 2010. 

[28] A. Aldroubi, A. G. Baskakov, and I. A. Krishtal, “Slanted matrices, Banach frames, and sampling,” J. Funct. Anal, vol. 255, pp. 1667-1691, 2008. 

[29] E. Cordero, K. Grochenig, F. Nicola, and L. Rodino, “Wiener algebras of Fourier integral operators,” J. Math. Pures AppL, vol. 99, pp. 219-233, 2013. 

[30] P. L. S0ndergaard, “Gabor frames by sampling and periodization,” Adv. Comput. Math., vol. 27, pp. 355-373, 2007. 

[31] A. J. E. M. Janssen, “The Zak transform: a signal transform for sampled time-continuous signals,” Philips J. Res., vol. 43, pp. 23-69, 1988. 

[32] C. deBoor, A Practical Guide to Splines. Springer, 1978. 

[33] O. Christensen and P. Massopust, “Exponential B-splines and the partition of unity property,” Adv. Comput. Math., vol. 37, pp. 301-318, 2012. 

[34] N. Dyn and A. Ron, “Recurrence relations for Tchebycheffian B-splines,” J. Anal. Math., vol. 51, pp. 118-135, 1988. 

[35] A. Ron, “Exponential box splines,” Constr. Approx., vol. 4, pp. 357-378, 1988. 

[36] T. Kloos, “Zeros of the Zak transform of totally positive functions,” J. Fourier Anal. AppL, DOI 10.1007/s00041-015-9402-5. 

[37] H. Bolcskei, F. Hlawatsch, “Discrete Zak transforms, polyphase transforms, and applications,” IEEE Trans. Sign. Proc., vol. 45, pp. 851-866, 1997. 

[38] S. Jaffard, “Proprietes des matrices ”bien localisees” pres de leur diagonale et quelques applications,” Ann. Inst. H. Poincare Anal. Non Lineaire, vol. 7, 
pp. 461^76, 1990. 

[39] A. G. Baskakov, “Wiener’s theorem and asymptotic estimates for elements of inverse matrices,” Funktsional. Anal, i Prilozhen, vol. 24, pp. 64-65, 1990. 

[40] R. Hagen, S. Roch, and B. Silbermann, C*-Algebras and Numerical Analysis. Monogr. Textbooks Pure AppL Math., vol. 236, Verlag C. H. Beck, 2001. 



